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Abstract. We consider a density estimation problem arising in nuclear physics. Gamma 
photons are impinging on a semiconductor detector, producing pulses of current. The 
integral of this pulse is equal to the total amount of charge created by the photon in 
the detector, which is linearly related to the photon energy. Because the inter-arrival of 
photons can be shorter than the charge collection time, pulses corresponding to different 
photons may overlap leading to a phenomenon known as pileup. The distortions on 
the photon energy spectrum estimate due to pileup become worse when the photon rate 
increases, making pileup correction techniques a must for high counting rate experiments. 

In this paper, we present a novel technique to correct pileup, which extends a method 
introduced in iHall and ParkI ll2004ll for the estimation of the service time from the busy 
period in M/G/oo models. It is based on a novel formula linking the joint distribution 
of the energy and duration of the cluster of pulses and the distribution of the energy of 
the photons. We then assess the performance of this estimator by providing an expression 
of its integrated square error. A Monte-Garlo experiment is presented to illustrate on 
practical examples the benefits of the pileup correction. 

Keywords: indirect observations; marked Poisson processes; nonlinear inverse problems; 
nonparametric density estimation 


1. Introduction 


We consider a problem occurring in nuclear spectroscopy. A radioactive source (a mix¬ 
ture of radionuclides) emits photons which impinge on a semiconductor detector. Photons 
(X and gamma rays) interact with the semiconductor crystal to produce electron-hole 
pairs. The migration of these pairs in the semiconductor produce a finite duration pulse 
of current. Under appropriate experimental conditions (ultra-pure crystal, low tempera¬ 
ture), the integral over time of this pulse of current corresponds to the total amount of 
electron-hole pairs created in the det ector, which is p roportiona l to the energy deposited 
in the semiconductor (see for instance iKno ll (|l989ll or Iheol (|l994l h In most classical semi¬ 
conductor radiation detectors, the pulse amplitudes are recorded and sorted to produce an 
histogram which is used as an estimate of the photon energy distribution (referred to in 
nuclear physics literature as energy spectrum). 
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(see Figure Q). ln~ga ):pnia ray specfr ometiy, Ttjiis phenomenon isj eterred to~as pileup. The 
pileup phenomenon induces a distortion of ihe acquired energy spectrum which becomes 
more severe as the incoming counting rate increases. i This problem li as been e xtensi vely 
studied in the held of nuclear instrumentation since the 1960’s (see Bristowl ( 199flll for 


a detailed revi ew of these early contributions; classical pileup rejection techniques are 
detailed in the ANSI nom] (199^). 
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Figure 1 . Illustration of the pile-up phenomenon; input signal with arrival 
times Tj, lengths Xj and energies Yj, j = n,... ,n + 2. Here X'^ = Xn, 
'^n ~ '^n, A^n+l = Tn+2 ~ Tn+l + Xn+2 and = Yn+l -|- Yn+2- 


In mathematical terms, the problem can be formalized as follows. Denote by {T^, k > 1} 
the sequence of arrival times of the photons, assumed to be the ordered points of an 
homogeneous Poisson process. The current intensity as a function of time can be modeled 
as a shot-noise process 

fF(f) = (1) 

k>l 

where {Fk{s), k > 1} are the contributions of each individual photon to the overall in¬ 
tensity. By analogy with queuing models, we call {W{t), f > 0} the workload process. 
The current pulses {Fk{s), k > 1} are assumed to be independent copies of a continuous 
time stochastic process {F{s), s > 0}. The pulse duration (the duration of the charge 

collection), dehned as X sup{f : F{t) > 0} is assumed to be hnite a.s. and the support 
of the path of F is assumed to be of the form [0,X] a.s., so that a busy period arrival 
corresponds to a pulse arrival and a pulse cannot belong to several busy periods. The 
integral of the pulse Y F{u) du is equal to the total amount of charge collected for a 

single photon. Under appropriate experimental condition, this quantity may be shown to 
be linearly related to the photon energy; for convenience Y is referred to as the energy in 
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the following. For the /c-th photon, we dehne the couple (Xfc, Yk) accordingly with respect 
to Ffe. The restriction of the workload process to a maximal segment where it is positive 
is referred to as a busy period, and where it is 0 as idle. In the coverage process literature, 
these quantities are also referred to as spacings and clumps. An idle period followed by a 
busy period is called a cycle. 

In our experimental setting, the sequence of pulse duration and energy {(X^, Y^), ^ > 1} 
is not directly observed. Instead, the only available data are the durations of the busy and 
idle periods and the total amounts of charge collected on busy periods. Dehne the on-off 
process 

Fi = , (2) 

k>l 


where {T^, A; > 1} is the ordered sequence of busy periods arrivals and {X'^, /c > 1} the 
corresponding sequence of durations. We further dehne, for all /c > 1, Yl ^(^) 

the total amount of charge of the /c-th busy period. Finally we denote by Zk the duration 
of the fc-th idle period, Zi = Ti and, for all k > 2, Zk = Tj. — (T^_i -|- X'f._i). We 
consider the problem of estimating the distribution of the photon energ y Y w i th n cycles 
{{Zk, X{, Yl), k = 1,... ,n} observed. In the terminology introduced bv lPvkel ( 1958h . this 
corresponds to a type II counter. 

The problem shares some similarity with service time distribution from busy and idle 
measurements in a M/G/oo model (see for instance Baccelli and Bremandl ( 2002(1 1. Note 
indeed that the M/G/oo model is a particular instance of the above setting, as it cor¬ 
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the cumulative dist r ibutio ns functions (cdfs) of the service time X and busy period X'. 
Bingham and Pitt3 ( 1999(1 derived from this formula an estimator of the service time dis¬ 


_ Y. Th ere exists a vast literature for this particu- 

see also nil (HHi)) has derived a closed-form relation linking 


tribution X, which they apply t o the study of biologi cal signals. An alternative estimator 
has been recently introduced in (Hall and ParlJ (jio^, in which a kernel-estimator of the 
probability density function (pdf) of X is derived in a nonparametric framework, together 
with a bound of the pointwise error. 

Although our estimator can be ap plied to the M/G/cxo framework (thus allowing a 
comparison with (Hall and ParlJ ((2004(1 in this special case), we stress the fact that we are 
dealing here simultaneously with durations and energies, without assuming any particular 
dependence structure between them. Secondly, the main emphasis in the photon problem 
consists in estimating the distribution of the photon energy and not the distribution of the 
duration, in sharp contrast with the M/G/cxd problem. 

The paper is organized as follows. We give the notations and main assumptions in Section 
(K and list the basic properties of the model. In Section (n( we present an inversion formula 
relating the Laplace transform of the cluster duration/energy to the Laplace transform of 
the density function of interest. We also derive an estimator of this function, which is based 
on an empirical version of the inversion formula and kernel smoothing. Our main result 
is presented in Section 01 showing that this estimator achieves standard minimax rates in 
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the sense of the Integrated Squared Error when the pulse duration is almost-surely upper 
bounded. The study of this error is detailed is Sectional Some applications and examples 
are shown in Section El Since the present paper is directed towards establishing a theory, 
pr actical aspects are no t discussed in much detail in the present contribution and we refer 
to iTrigano et ah ( 200511 for a thorough discussion of the implementation and applications 
to real data. Proofs of the different propositions are presented in appendix. 


2. Notations and main assumptions 

All along the paper, we suppose that 

(H-1) {Tfc, A; > 1} is the ordered sequence of the points of a homogenous Poisson process 
on the positive half-line with intensity A. 

(H-2) {(X, Y), (Xfc, Yfc), /c > 1} is a sequence of independent and identically distributed 
(0, cx))^-valued random variables with probability distribution denoted by P and 
independent of {T^, /c > 1}. In addition, E[X] and E[y] are hnite. 

In other words, {(Tfc,Xfc,Yfc), A; > 1} is a Poisson point process with control measure 
ALeb (8) P, where Leb denotes the Lebesgue measure on the positive half-line. Let us recall 
a few basic properties satished under this assumption by the sequence {(2’^, X{,, A: > 1} 
dehned in the introduction. By the lack of memory property of the exponential distribution, 
the idle periods are independent and identically distributed with common exponential 
distribution with parameter A. Moreover they are independent of the busy periods, which 
also are independent and identically distributedWe denote by (X', Y') a couple having 
the same distribution as the variables of the sequence {(X^,^^'), A: > 1} and by P' its 
probability measure. Using that E[X] and E[y] are hnite, it is easily shown that 

E[X'] = {exp(AE[X]) - 1}/A 
E[r'] = E[y] exp(AE[X]) . 

Our goal is the nonparametric estimation of the distribution of Y ; hence we assume that 

(H-3) Y admits a probability density function denoted by m, i.e. P{dx, dy) = 

m{y)\jeh{dy). 

As mentioned in Section [H the marks {{X^^Yk), A; > 1} are not directly observed but, 
instead, we observe the sequence {(T^, X[, X^'), A; = 1,... ,?7,}, i.e. the arrival times, dura¬ 
tion and integrated energy of the successive busy periods. These quantities are recursively 
dehned as follows. Let T[ = Ti and for all A; > 2, 


T> 

-‘-k 


inf iTr.T,> 





(3) 
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for all k > 1, 


^k = ^ max {Ti + XJ - , 

i>l 


( 4 ) 


Remark 2.1. In this paper, it is assumed that the experiment consists in collecting a number 
n of cycles. Hence, the total duration of the experiment is equal to T^ + X^ and is therefore 
random. A classical renewal argument shows that, as n —> oo, + X^jn converges a.s. 
to the mean duration of a cycle, 1/A + E[X] = Aexp(AE[X]). Another approach, which is 
more sensible in certain scenarios, is to consider that the total duration of the experiment 
is given, say equal to T. In this case, the number of cycles is random, equal to the renewal 
process of the busy cycles, Xj = '^T=i ^{'^k — T}. As T oc, the Blackwell theorem 
shows that Nj/T —> 1/A exp(AE[X]), showing that the asymptotic theory in both cases 
can be easily related. 


3. Inversion formula and estimation 


Let P be a probability measure on M x M equipped with the Borel a-algebra; for all 
(s,p) G X C"*", where C"*" = {z E C, Re(. 2 ) > 0}, we dehne its Laplace transform (or 
moment generating function) CP as: 


CP{s,p) 



P{du, dv) . 


The following theorem provides a relation between the joint distribution of the individual 
pulses energies and durations P and the moment-generating function of the distribution of 
the energies and durations of the busy periods CP'] this key relation will be used to derive 
an estimator of m. 


Theorem 3.1. Under Assumptions\(H-T][\(H-2)\ for all {s,p) G C"*" x 


/•+00 

ir / N nj ACP'(s,p) 1 

{a(u,p) l}du = 

s -|- A s -|- A — XCP'{s, p) 

(5) 

J u=0 

where 

a{u,p) exp (AE[e“^^(M — X)+]) . 

(6) 

Proof. See Section iDl 


□ 

Remark 3.1. Observe that the integral in (jnj can be replaced by ^ since. 

in dn}. 


a{u,p) = 0 for M < 0. Moreover, from Q, we trivially get \a{u,p)\ < exp(Au) for Re(p) > 0; 
hence this integral is well defined for Re(s) > 0 and Re(p) > 0. 

The relation © is rather involves and it is perhaps not immediately obvious to see how 
this relation may yield to an estimator of the distribution of the energy. By logarithmic 
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differentiation with respect to x, implies 

loga(a;,p) = AE[e“^^l(X < x)] . (7) 

ox 

We consider a kernel function K that integrates to 1 and denote by K* its Fourier 
transform, K*{y) = dy, so that K*{0) = 1. We further assume that K* is 

integrable, so that, for any y eM., 

K{y) = ^ r K^vy'^ydu. 


Hence, from dlD and Fubini’s theorem, we have, for any bandwidth parameter h > 0 and 
all 1 / G M, 


1 


1 d 


log a{x, iu)K* {hu)e^'^ydu 


= E 

= E 


A dx 

— / K*{hu)e^‘'^y-^H{X < x)du 

27r J _r^ 


K 


y-Y 

h 


1{X < x) 


( 8 ) 

(9) 


Taking the limits x —>■ oo and h —> 0 in the previous equation leads to the following explicit 
inversion formula which will be used to derive our estimator. For any continuity point y 
of the density m, we have 

f 1 7 +“ 1 <9 1 

m(y) = lim lim <— / ——log a(x, ii^) K*(hiy) diy > . (10) 

A ox J 


We now observe that for any p G C'*', the RHS of © is integrable on a line {c + io;, a; G M} 
where c is an arbitrary positive number. By inverting the Laplace transform, implies 
that, for all p G C"*" and x G M+, 


a{x,p) = 1 + 

^ 7+“ CP'{c + iuj,p) _ 

27r Jc T luj T A c T iu T A — {c T ico’, p) 

Our estimator of m is based on m and m but we need first to estimate A, the intensity 
of the underlying Poisson process. Since the idle periods are independent and identically 
distributed according to an exponential distribution with intensity A, we use maximum- 
likelihood estimator based on the durations of the idle periods {Zk, k = 1,..., n}, namely. 


• ( 12 ) 

The function a(x, w) can be estimated from {(X^, Y)f), = 1,..., n} by plugging in (fTT| an 

estimate of the Laplace transform CP' of the joint distribution of the busy period duration 
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and energy. More precisely, let P'n be the associated empirical measure: for any bivariate 
measurable function g, we denote by 

P'ng= [jg{x,y)P'n{dx,dy) = ^'^g{Xk,Yk) . 

^ k=i 

We consider the following estimator 


an{x, w) = 1+ 


^+oo 




XnCP'n{c + icn, iu) 


3(A+c+itij)a: 


c + icn + A„ c + io; + An — CP'n{c + icn, iu) 


where 


CP'nic + icn, ii^) CP'nic + ictJ, ii^) = — e 


-{c+iw)X'^-ivYl 


duj . (13) 


(14) 


k=l 


In practice, the numerical computation of this integral (and also the one in (ETll below) 
can be done by using efficient numerical packages (we can refer to iGantschil (|l997ll for an 
overview of numerical integration methods). Since the integrand is inhnitely differentiable 
and has a modulus decaying as when uj —> ±cx), the errors in computing this integrals 

numerically can be made arbitrary small. The numerical error will thus not be taken into 
account here for brevity. 

In order to estimate \~^d\oga/dx, we also need to estimate the partial derivative da/dx. 
Because the function x ^ a{x, iv) (see (HU) is dehned as an inverse Fourier transform of an 
integrable function, it is tempting to estimate its partial derivative simply by multiplying 
by a factor A + c + itn its Fourier transform prior to inversion. This approach however 
is not directly applicable, because multiplying the integrand by a; in (FTTl) leads to a non 
absolutely convergent integral. As observed bv iHall and Parl3 (200^ in a related problem, 
it is possible to get rid of this difficulty by hnding an explicit expression of the singular 
part of this function, which can be computed and estimated. Note hrst that, for any s 
and p with non-negative real parts, \CP\s,p)\ < 1; on the other hand, Re(s) > 0 implies 
|A/(s -I- A)| < 1. Therefore, for all (a;, z/) e M x M, 


E( 


c “h i(jj -f" A — Xj(1P'(c iuj^ iz/) c itu -|- A ^ V 

^ ' n>0 

Using the latter equation, we obtain 

XCP'{c + iu;/w) 1 

c “h iuj H" A c P iuj “h A — XHP^i^c H- iuj^ iz^) 


XCP\c + iu;, iu] 
A c “t“ ico^ 


= x4i(a;, iz^) + ^ 2 ( 0 ;, iz/) (15) 















E. MOULINES, F. ROUEFF, A. SOULOUMIAC AND T. TRIGANO 


where we have dehned 

A2{uj,w) 


XCP\c + io;, w) 

(c + iu; + A)2 

{A£P'(c + iu;, 1 

(c "h iu -|- A)^ c iu A" A — XPP'{c P icu, iz/) 


It is easily seen that the functions uj i—*• Ak{oJ,iv), k = 1,2 are integrable. Hence we may 
dehne, for k = 1,2, and all real numbers x and u. 


/ . s def 1 

ak{x,iu) = 


2TI\ J 


Ak{u,-w) eA+c+iPx 


(16) 


and therefore, using m and m, a^x,!!^) = 1 + Aai(a:, ii/) + Aa 2 (x, ii/), which hnally yields 


log a{x,w) = ——— 
A ox a{x,w) 


dai da2 
dx dx 


X, w 


(17) 


Recall that the moment generating function of a gamma distribution with shape parameter 
2 and scale parameter A is given by a: i—*• A^/(A — x)'^. It follows that, for all m G M, 


„(c+itj)u 


du = u+e 


—Xu 


‘2tt Ju}=-oo (A + c + iujy 
Using Fubini’s theorem and this equation, we obtain, for all real numbers x and v. 


ai{x,IV) = 

= e^^E 


]gJe-((c+it^)A'+ii/y')jg(A+c+itj)a 


27r / 

J uj = — oo 


= E 


(A + c + iujY 

e-'^’^' {x - 

>\ \X’-iuY 


doj 


{x-X')+e^ 


and, differentiating this latter expression w.r.t. x, we obtain 


dai 

dx 


{x, iv) = E 1{X' < x) e 


\X'-ivY' 


( 18 ) 


On the other hand, note that 1 ^ 2 ( 0 ;, \v) \ = 0(|a;| ^) as a; —*■ ±cx), the derivative of 02 can 
(and will) be computed by multiplying the integrand in (fTKI) by A + c + iu;, namely. 


9(^2 / ■ \ ^ 


^ + CXD 


X 


{CP'{c + iu;, iz^)}^ 


^(A+c+ia;)fc 


c -f- luj ~h A c -|“ luj -f- A — A>C-P^(c -|- ico^, iz/) 


duj . (19) 
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Eq. and m then yield the following estimators for dak/dx, k = 1,2, 


h,n{x, = - y] l(y < x) 

^ T p{c+\n)x 

hAx,^) = -;-X 




( 20 ) 


2ir 


^H-oo 


{CP'n{c + lOJ,lv)Y 


1-00 c + io; + c + icn + An — An£P'n(c + icn, ii^) 


duj 


( 21 ) 


where An and CP'n are given respectively by m and (HI. From (HUD and (HID, we hnally 
dehne the following estimator for the energy distribution density function: 




271 


^+oo 

h,n + l2,n / . N 

/ 

[X, IP) 

/ —oo 

dn 


K\hv)e^'^y dv , 


( 22 ) 


where dn, h,n and l2,n are respectively dehned in HSD, (EHD and m- 


4. Main result 


We denote respectively by || ■ ||oo, || ■ Ih and || ■ ||w(/ 3 ) fhe inhnite norm, the L^-norm and 
the Sobolev norm of exponent (d, that is the norm endowing the Sobolev space 

>V(^) = G L2(M) ; = J (1 + \iy\f^\g*{iy)\'^du < ooj , 

where g* denotes the Fourier transform of g. Consider the following assumption on the 
kernel. 


(H-4) K* has a compact support, and there exists constants Ck > 0 and I > {3. such 
that for all i/ G M, 


\1-K*{u)\<Ck 



\u\ 

1 

(1 + 

\u\y 


We may now state the main result of this section, which establish the rate of convergence 
of the integrated square error. 

Theorem 4.1. Let (3, C and x be positive numbers. Assume w^TWim and suppose 
that X < X a.s. and ||m||>v(^) < C. Then there exists C' > 0 only depending on K, X, c, 
(3, X and C such that, for all M > 0, 

limsupP(n^/(^+2/3 )||> M) < C M-^ . (23) 

n—*^oo 

where hn = ^ 

Proof. See Section [HI □ 
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Remark 4.1. In the application we have considered, the condition X < x assumption is 
always satished. Indeed, the pulse duration corresponds to the duration of the charge 
collection, and therefore to the lifetime of the pairs of electron-holes in the semiconductor 
detector. This lifetime is always hnite and depends primarily of the geometry of the 
detector. 


Remark 4.2. However, the condition X < x a.s. can actually be circumvented if, at hxed 
X > 0, one considers mx,hn,n as an estimator of dehned as the density of the measure 
/JLowhich is always dehned under Assumption |(H-3)j 

Remark 4.3. If X and Y are independent, then mx{y) = m{y)F{X < x) so that, for all x 
such that P(X < x) > 0, we obtain an estimator of m up to a multiplicative constant. 

Remark 4.4. In the M/G/oo case, i.e. ii X = Y a.s., = m l[o,x]- Hence, since 


'^x,h„,n ^)l[0,x]||2 — ^ l[0,x]) 


|2 

I2 ’ 


our results apply to the locally integrated error for estimating m l[o,x]- As a comparison, 
the rate of our estim ator is given by the s moothness of whereas the rate of the 

estimator proposed in iHall and Parkf (|2nn4 ) for estimating the time service density is given 
by the smoothness of the pdf of X' (see ( Hall and Parkl. 2004, Eq. (3.7))). 


Remark 4.5. The estimators in (ESI) are functions of {(Z^, Y)'), /c = 1, ..., n}, where n 

is the number of observed cycles. For t G M+, denote by A/) the renewal process associated 

to the arrivals of the photons. A/) '= ^{'^k < t}. The number of arrivals during 

n cycles is equal to h = and is therefore random. As n tends to inhnity, (T,( -|- 

X'^/n converges almost surely to the mean of the cycle duration, exp(AE[X])/A and it 
can be easily shown that the u-th return to an idle period (that is, + X'^) is a stopping 
time with respect to the natural history of A/). Therefore, by the Blackwell theorem, 
A/r^+x4/(7(( -|- X'^) converges to A. Therefore, n/n = AfT^+x'^ln converges almost surely 
to exp(AE[X]). It is well known that the minimax integrated rate for estimating m from 
{Yfc, /c = 1,... ,h} with m in a /9-Sobolev ball is the only non-standard feature 

being that the density estimator is calculated by using a random number of data, which 
does not alter the density’s estimator hrst-order property. Since h/n converges almost 
surely to a constant. Theorem 14.II shows that the rates achieved by our estimator is the 
minimax integrated rate. 


5. Decomposition of the error 

We give in this section theoretical results for the proposed estimators. We hrst introduce 
auxiliary variables, which will be used in the proof of the main theorem. For any positive 
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numbers W, x and A, define 

An(Vb) = sup 

{uj,u)e[-w,w]‘2 


En{W-,x, A) 


def 

= sup 
ue[-w,w] 


\CP'{c + icT, iv) — CP'n{c + icT, iu) I ; 

1[0,x](m) e-‘""(P' - T'n)idu, dy) 


Proposition 15.II provides bounds for the random variables A„ and 


Proposition 5.1. Assume yjnwim Then Mi = E(max{X', Y'}) is finite and the 
following inequalities hold for all e > 0, r > 0 and IP > 1; 


^ 4 r ilP 1 

P(|A„(M')|>£)<^-!.+ 



(24) 


^ ~ At Ml ( IP A ( n£^ \ 

sup n\En{W-,x,X)\ >e)< -1 + 1 + ^ exp . (25) 

a.,A>0 ^ V ^ / V 16 y 

Proof See Appendix O □ 

Since our estimate depends on A^ and CP'n, we introduce auxiliary functions to ex¬ 
hibit both dependencies. Define the following functions depending on h, x, A and on any 
probability measure P : 

~ ~ Hef e^'^+^A r+°° 

d{x,ih';X,P) = 1-1-/ <h(c-|-iu;, ip ; X, CP) duj (26) 

2vr J_oo 


/i(x,ip; A, P) 


del 


l{u<.}e^"-"'^P{du,dv) 


J2(x, P;A, P) = 


def 6 


(A+c)x p+oo 


271 


and define 


m(y; x, h, X, P) = ^ 


CP{c + iu;, iv) <h(c + iu;, ii^ ; A, CP) du 


K*{hiy)fi'^y du (27) 


r+oc 

Il+ I 2 ( . T 5\ 

„ (x, IZ/; A, P) 

J —oo 

a 


whenever the integral is well defined. Hence, by (ITH|l . (On|l and (j^H), for 

*=1,2, 


d{x, w] A, P') = a{x, xp) 

and 


doi 

(28) 

d{x, iv] Xn, P'n) = dn{x, iv) 

and 


= li^nixfiv) 

(29) 
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bi{,y) ‘^= m{y) -E 


and m^xniy) = x, h, Xn, P'n)- Now define 

b- 

y-Y 


62 ( 1 /) =E 


K {y) = miy; x, h, A, P') - rh{y, x, h, A„, P') ; 


1 


y-Y 

h 


-m{y,x,h,X,P') ; 


V 2 {y) rh{y, x, h, Xn, P') - mx,h,n{y) 


so that, by definition, 


(30) 

(31) 

(32) 

(33) 


'IT^x,h,n — m = &1 + 62 + hi + V 2 . (34) 

In this decomposition, hi and 62 are deterministic functions and Vi, V 2 are random pro¬ 
cesses. We now provide bounds for these quantities in the L? sense. 


Theorem 5.1. Let jS, x and h be positive numbers and n be a positive integer. 
sume WTwim If m & hV(/3), then we have 

||^l||2 — ll"^llw(/3) i (35) 

\\b 2 \\l<\\K\\lh-^¥[X>x\. (36) 

Moreover, there exist positive constants M and y only depending on c and X such that the 
two following assertions hold. 

(i) We have 

III/ 1 II 2 < M‘^\\K\\l{l + xfh-^e^^^^^^^^(Xn-Xf (37) 


= <1 |A„ - A| < r/ (1 + a;)-i 


on the event 


(a) For allW >1 such that [—Wh,Wh] contains the support of K*, we have 


11^2112 < \\K\\l h-^ \An{W) + W + E„(W; x, A^ 


12 


(38) 

(39) 


on the event Ei intersected with the event 

E2 |a„(W) + W-^ < . 

Proof. See Appendix 1X1 


(40) 

□ 


In this result, hi is the usual bias in kernel nonparametric estimation; 62 is a non-usual 
bias term which only vanishes when X is bounded, it correspond to the fact that the 
limit a: —> cx) is not attained in (Uni); the fluctuation term Vi accounts for the error in the 
estimation of A by Xn and is of the order h~^y/n for hxed x and V 2 accounts for the error 
in the estimation of CP' by CP'n and, by using Proposition 15.11 it can be shown to be 
“almost” of the order h~^y/n for W chosen to diverge quickly enough with respect to n. 
The events Ei and E 2 have probability tending to 1 as n tend to inhnity; they are induced 
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by the fraction present in the definition (|22D of the estimator as they primarily avoid the 
denominator approaching zero. 

We now give a result on the consistency of our estimator, and also on a rate of conver¬ 
gence, based on Theorem 15.11 and Proposition 15.11 bv imposing a superexponential tail for 

X. 

Corollary 5.1. Let (3 > 0 and 7 > 1. Assume and suppose that m G W(/9) 

and F[X > x] = 0(e“l^l^). Then, for all e > 0, as n ^ -|-oo, 

= , (41) 

where hn x and Xn x (logn)'^' with 7 ' G ( 7 ”^, 1 ). 

Proof. We set Wn "= n. By Proposition 15.11 we get by choosing e = C(log(n)/n)^/^ and 
r = : 

P(|A„(1W)| > C(log(n)/n)^/^) < ^log"^/^(n) + (1 + 

which tends to 0 as n —> cxo for C > 32. Hence, 

|A„(HA)| = Op{(log(n)/n)^/^} . 

Similarly, since is independent of k = 1 ,...,? 7 ,}, we get x„, A„)| = 

Op{log(n)/n)^A Observe that, for any 5i > 0, (52 > 0 and e > 0, exp(( 52 X„) = o(n'^). 

Since A„ = A -|- Op(n“^/^) and |A„(HA) -|- Wp^\ = Op((log(n)/n)^/^), Ei and E 2 have a 
probability tending to one, so that the bounds of Theorem 15.11 finally gives, for any e > 0, 

m\=0^{{hnny-^/^) , * = 1 , 2 . 

Now using the superexponential tail assumption for X, we have P(A > Xn) = 0(exp{— log"’'"’' (**)}) 
o{ny for all e > 0 and the result follows. □ 

As seen from m, the estimator almost achieves the standard nonparametric minimax 
rate that one would obtain by observing {(A^, W), k = 1,... ,n} directly. If X 

is bounded, then the rate can be made more precise as in Theorem 14.11 by taking x equal 
to an upper bound for X (so that 62 = 0) and hn x one easily gets from the 

above proof that 

\\m - = Of (log(n)n"2^/(^+^^)) , (42) 

thus a lost of log(n) in comparison with the claimed rate. This log(n) can in fact be 
removed as shown in Appendix IHl 


6. Applications — Discussion 


The present paper is directed towards the construction of an estimator and deriving 
elements of its asymptotic theory. We will therefore sati sfy ourselves by providing simple 
examples and will refer the reader to TriEfano et ah IJiooi) for an in-depth discussion of the 
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selection of the setting parameters {e.g. the kernel bandwidth, the truncation bound, etc) 
and the analysis of many different data sets. 

We hrst consider a simple simulated data set. Samples are drawn according to the 
bimodal density 

fi^yV) = A/ 20 , 3 (^) X (0.6A/ioo,6(|/) + 0.4A/i3o,9(|/)) , (43) 

where denotes the gaussian distribution of mean a and standard deviation b truncated 
to R+; The intensity of the Poisson process is set to A = 0.04. Figure [21-(a) shows the 
true density and a kernel estimate of the marginal of the pileup distribution, based on 
10® samples. Figure EKb) displays the difference between the true and the estimated 
density, obtained using the kernel bandwidth h = 2 and the upper bound x = 80. We 
see that the estimated energy distribution mx,h,n captures most of the important features 
of the original distribution m. Note that the second mode of the original density is well 
recovered after the pileup correction, whereas it is severely distorted in the absence of 
any processing. The fake modes that appear at energies 200 and 230 are totally removed. 
This is nevertheless a toy example, since we pointed out that in our application X and Y 



(a) (b) 


Figure 2. (a): marginal density of the energy for the pdf (dashed) 
and kernel density estimator of the pileup observations based on 10® samples 
(solid), (b): marginal density of the energy for the pdf (dashed) and 
estimator rhx,h,n (solid). 


were not independent. Numerical values of the mean integrated squared error (MISE) are 
presented in Tabled for a fixed bandwidth parameter h = 2.0 and different values of n, c 
and X. It shows that c has little influence on the error. This is hardly surprising, since the 
Bromwich integral used to comp ute the inverse Laplace transform does not theoretically 
depend on the choice of c (see e.g. Doetschl (19^). Concerning the influence of x, knowing 
that X has distribution A/ 20 , 3 , “reasonable” values (displayed in the three hrst rows) all 
give equivalently good results but the last row shows that the “naive” data-driven choice 
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X = maxj<„ X- significantly deteriorates the estimate. Indeed, in view of the upper bounds 
of Theorem o on the one hand, choosing x too large does not ensure that the variance 
term Vi and V 2 are controlled, since in this case the conditions and (P|) may not be 
satished; on the other hand, x too small introduces a bias in () 22 | 1 . since the control of the 
bias term 62 is not guaranteed in that case. 



Table 1. Mean Integrated Square Error Monte-Carlo estimates as n, c or x varies. 


We now present some results using a more realistic model of the energy distribution of 
the Cesium 137 radionuclide (including Compton effect). We draw n = 500000 samples of 
{X, Y) using the adaptive rejection sampling algorithm, according to the following density: 


f{x,y) = m{y) x fx\Y{x\y) , 

where m is represented by the dotted plot of Figure 01 (a) and the conditional distribution 
fx\Y{'\y) is a Gamma distribution with unit scale parameter, shape parameter equal to 
2 -|- I//1024 and truncated at T = 4 + 1//2048; the number of samples may appear to be 
large, but such large number are commonly used in nuclear spectrometry, especially when 
active sources are measured. Figure 01 (a) also shows the pileup distribution (solid curve), 
based on the observations Yjf, k = 1 ,... ,n to illustrate the difference with m; note that 
the Compton continuum (which is the smooth part of the density on the left of the spike) 
is also distorted, since electrical pulses generated by Compton photons are also susceptible 
to overlap. Figure El (b) illustrates the behavior of our estimator. We observe that the 
pileup effect is well corrected. 

We now briefly discuss on the choice of the bandwidth parameter h. In standard nonpara- 
metric estimation, there are several data-driven ways of choosing a bandwidth parameter. 
It is not yet clear how these methods can be adapted to this non-standard density esti¬ 
mation scenario, except in special cases. For instance, a possible approach would then 
consist in using an automatic bandwidth selector (such as cross validation) on the observa¬ 
tions {Yl, k = 1,..., n}, and use the obtained optimal bandwidth for the estimator mx,h,n- 
Further insights on the data-driven choices of c, x and h and d i scussio n of the practical 
applications can be found in the companion paper iTrigano et al. f 20051) . 


Appendix A. Proof of Theorem 15.11 
The following lemma will be used repeatedly : 
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Figure 3. Energy spectrum of the Cs 137 element — (a) Ideal probability 
density function (dotted curve) and kernel estimates of the pileup distribu¬ 
tion (plain curve); (b) Estimate rhx^h,n 


Lemma A.l. Let c > 0 and rjo > 0. for any complex valued functions Zi and Z 2 satisfying 

sup \zi{ijj)\<l , (44) 

a;GM,2=l,2 

let z = (zi, Z 2 ) and denote by the function defined on M+ x M. by 

^ ^i(^) 

^^- . 

(c -|- ia; -|- A)(c -|- iu; -|- A — \z 2 {uj)) 

Then the following assertions hold : 

(i) The function A 1 —*• 4/^(A, a;) do; is continuously differentiable on M_|_ and its 

derivative is bounded independently of z over A G [ 0 , 170 ]- 

(ii) There exists K > Q only depending on c and rjo such as, for any PF > 1 or IE = cx) 
and any function z = (ii,i 2 ) also satisfying 


sup 

Ae[o,T;o] 


'•+00 


4/^(A, uj) duj 


b+oo 


\l/i(A, uj) duj 


< K 


( max sup \zi{uj) 

*=^’2 0J&[-W,W] 
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Proof. For all ci; in M and A G [0, r^o], by using (jUJ, we obtain 


a3,d'(A,a;) 


< 


3(c+ |a;|) + 47^0 


2 • 


(45) 


(c^ + a;2) (c + A)2 + a;^ — A 

From Jensen inequality, one can easily show that for all a in [0,1] and A G [0,77o]) 

\J {c + A)^ + o;^ > (c + X)y/a +\u:\\/l — a . (46) 

Choosing a close enough to 1 so that (-y/a — l)?7o + c-y/a > 0, (PH|) and yields to 

3(c + |a;|) + dr^o 




< 


(c^ + a;2) [{^/a - l)r]o + c^/a + \uj\^/l - a) 


2 ’ 


which is valid independently of A in [0, rjo], and whose RHS is integrable over uj in M, hence 
(i). For showing (ii), observe that for all A in [0,?7o] 


^,(A,o;) - 4^5 (A,cu) = 


(2:1(0;) - zi{uj)) 


(c J- io; J- A — A2:2 (o;))(c J- io; J- A — Az2(o;)) 

A ( Ai ( a ;) 2 ; 2 ( a ;) - Z2{uj)zi{uj)) 


+ 


(c J- io; J- A)(c J- io; -|- A — Az2(o;))(c -|- io; R A — Az2(o;)) 

Using again iSD and since, using (jUJ, \z1Z2 — Z2Zi\ = \zi{z2 — Z2) — Z2(zi — zi)| < 
1^:2 — ^2! + l^i — Ji|, for a chosen as above and for each A in [ 0 ,r^o],, we have 


|4^,(A,o;) - 4^5 (A,o;)| < 


(c + |o;| + 3rjo) 


- o maxRho;) — Zi(u)\ . (47) 

\/c^ + o;2 ((i/a - l)?7o + Cy/a + |o;|\/l - a) *=^’2 

Observe that in the RHS above, the fraction is integrable over o; G M and is equivalent 
to [(1 — q:)|o;|]“^ as |o;| —> 00 . Consequently, there exists constants Ki and K 2 depending 
only on c and r/o such as 


r+W 

l-W 


- 4'i(A,o;)| duj 

< AJ max sup |2;i(o;) — Ai(o;)| 


U1&[-W,W] 

and, since maxj=i^2 \zi{uj) — Zi{uj)\ < 2, 




A/r^^(cv/a + {^/a - 1)770) 
4^,(A,o;)-^2(A,o;)|dco<iC2hF“^ , 


hence (ii). 


□ 
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Bound for hi. Observe that hi is the usual bias in nonparametric kernel estimation. The 
bound of the integrated error is classically given, for density in a Sobolev space, by 


II61II2 = / \l — K*{hv)\^ \m*(y)f‘ dv < \\rn\ 

J —00 

which shows (jsni). 

Bound for 62. By 0, |T7j) . (P7jl and (P|). we hnd 


W(/3) 


hiv) = E 




An application of the Cauchy-Schwarz Inequality yields 

Bound for Vi. We will show below that there exist positive constants M and rj such that, 
on El (as defined in dSHD), 


sup 

i/gR 


d 

'h + h' 

d\ 

d 


{x,iu;Xn,P') 




(48) 


Using (j77|) and (ES), the Parseval Theorem and the latter relation imply 

/ oo ^2 

M{l + x) - A) \K*{hu)\^ dv 

-OO 

on the event T'l, which yields (EZD- Hence it remains to show ()48j) . 

First observe that, by definition of /i, one gets trivially, for all A > 0, 


|Ji(a:, iz/; A, P')| < e"^"^ and 


A 


/i(x, iz/; A, P') < xe^^ 


Inserting into ((2H1), we get, for all z/ G M, 

|a(a:, iz/; A, P')| = exp (AE [cos(z/y) (x —X)+]) G [e“^'^,e^'^] 


(49) 


(50) 


Let r^o > 0 to be chosen later. From (I20|) . Lemma mu shows that, there exists a constant 
Ml only depending on A, c and rjo such that, for all A in [0, A + tiq] and z^ G M, 

{d{x, iz/; A, P') - l)e-(^+")^ - {d{x, iz/; A, P') - l)e-(^+")^ < Mi|A - A| . (51) 

From CT and (ED and since, for all real z/, |1 — e^| < |2/|el^(, we get, for all A in [0, A + z^o] 
and zz G M, 

\d\{x, izz; A, P') > \d\{x, izz; A, P') - 1| - Mi |A - A| , 

|a|(x,iz/;A,P0 >e(^-2^)^-[Mie(^+")^ + xe>^-^l^]|A-A| . 


hence 
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Note that, taking rjQ = c and M[ = Mi V 1, the term between brackets is at most 
e('^-2N*(l _|_ ^'^q(c+ 2 X)x ^ g gQ gg^^ Qj^ |;]^g gyent El with 

V < Vi = {Vo A (M0-V2}, 

|a|(a;,iz/;A„,P')>^e(^"-2")". (52) 

From (j2b|) and using similar bounds as in Lemma lA.ll one easily shows that, for some 
constant M2 only depending on A, c and rji, for all A G M such that |A — A| <171, 

\a{x, IV] A, P') I < M 2 

1/2(0:, iv] A, P')\ < M 2 and 

Collecting ii, (El and the two last displayed bounds shows that (HH|) holds on Ei, for 
any rj < rji. 

Bound for V2. Since the support of K* is included in [—Wh,Wh], By Parseval Theo¬ 
rem, (El and (El, the claimed bound is implied by 


\ TD'\ 

^(x, m;A,P) 
oX 

^(x, m;A,P) 
oX 


< M2 (1 + x) , (53) 

< M2 (1 + x) . (54) 


{x, \V] An, P') - iv] An, P'n) 


sup 

\u\<W ®. ® 

< M 

which we now show. Using El, we may write 

yg(c+A)x ^+00 


lXn{W)+W-^ +En{W]xAn) , (55) 


|a(o:, \v] A, P) — a{x, iv, A, P') \ = 


271 


[Ti(A,a;) - T^(A,u;)] du 


(56) 


where T is defined in Lemma EH and where the complex functions 5 and z are defined 


as 


z{uj) {P‘^^CP'{c + icT, iz/); CP'{c + iu, iv)) 


and 


z{uj) (e'‘^'^£P(c -I- io;, iz^); £P(c -|- icr, iv)) . 

Using El and assertion (ii) of LemmaEH there exists Mi > 0 such that, for all A < A-l-r^i, 


sup \a{x,iv; X, P') 

\u\<W 


d{x, iv, A, P'n) I < Ml e('=+")" ( An(lU) + — 


It is also clear that for all A < A -|- r^i, 

sup |Ji(x,iu; A, FV) - h(x,u^; A, F')| < e^En(W;X,x) 

\i>\<W 


(57) 


( 58 ) 
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and 

with 

and 


1/2(0:, iu; A, P) - hix, iu; A, P') \ = 


^(c+A)* r+00 


27r 


[d^2(A,a;) - dxj 


z{uj) + icn, CP'{c + icn, iv)) 


z{ijj) {P'^^{CP{c + icn, iz^))^; CP{c + io;, iz^)) . 
Consequently, using assertion (ii) of Lemma fA.il we have for all A < A + r^i, 

sup |J2(o:,iz/;A,:P„)-J2(o:,iz/;A,P')l (K{W) + ^r 

\v\<W \ rV 


( 59 ) 


We now derive a lower bound for a„(a:, iz/) = 0(0:, iz^; A„, P'n); By we get 

inf |a|(a:,iz/)> inf |a|(o:, m; A„, P') - ( A„(W) + 4 ; 

\u\<W \y\<W \ W 

Recall that Pi and P2 are dehned in and iOD respectively. Using (EH); which holds 
on Pi for any rj < rji, we get, on Pi fl P2, 


inf 

u\<W 



1 

n 1 

a(x, iz/; A„, P'„) a„(x, iz/) 

> - 
“ 2 

'to 1 1 
1 

03 

1 _ 


^(2U-4A)x 


( 60 ) 


Hence we set 77 ( 4 Mi) ^ A 771, so that the term between brackets is at least 1 / 4 . Finally, 

using that, for all complex number x, y, z, x', y', z', 

X + y x' + y' {z' — z){x + y) + z{x — x') + z{y — y') 


z 


z' 


zz' 


and collecting (jlHl), (ESI), (EH), (EH), (ESI), (ESI) (jSSl) leads to ESD- 


Appendix B. Proof of Theorem 14.11 

In this section we denote by 771, M* and Cj, 7 = 0,1, 2,... some positive constants only 
depending on ||m||>v(/3), P, A, c and x. We will also use the notations introduced Section 1X1 
As shown in this section, we have II/1II2 < and since P(P > x) = 0, we 

have 62 = 0. 

By ®, because E[X] < cx), it is easily seen that |a(a:, in)! > e This can be used in 
the ratio appearing in to lower bound a„ for n large by using that an{x, iz/) converges 
to a(x, iz/). However this will not allow bounds of the ratio in the mean square sense. For 
obtaining mean square error estimates, we consider the following modihed estimator which 
(artihcially) circumvent this difficulty. Let rjo > \ and denote by An the set 

An = {A„ < 770} n I inf \an\{x, iz/) > ^ exp(-Ana:) I , 

I Zxni'GSupp(A*) 5 J 
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where Supp(J’r*) denotes the (compact) support of K*. Dehne 

'^x,h,n(jj') ^An ■ 

We will show that 

sup - m\\l< Co . 

n>l 

Let V3 be the random process 

1/3(7/) = 7h{y] X, hn, A, P') - rn^,h„,n{y) > 

so that 

\\m^,h,n - m\\l < \\bi\\l + llmll^l^c + HI/sHsIa^ 

We will show that there exist Ci > 0 such that, for n large enough, 

nK)<Cin-'; 

E|IIU||p^J<Ci(/!„n)-‘. 


(61) 

(62) 


(63) 

(64) 


Since ||i>i||2 < W^Wwlp) ll^lle l|nt||w(/3)5 ihe three last displays yield the 

bound (El- The bound (El then follows by writing 


p(^/3/(l+2/3)||-^^^^^^ - m \\2 > M) 

< _ ^||2 > M) + P(Ai:) < CoM-^ + , 

where we applied the Markov Inequality, dSl and (El. It now remains to show El 
and (El- 

Proof of bound (j63jl . We set Wn *= n, so that for n large enough, hnV € Supp(iL*) 
implies |z/| < Wn- As in (El, we have, on Ei fl £'2, 

inf |a„|((z:,iz/) > _ 


Hence the intersection of {A„ < rjo}, Ei, £2 and {exp((A„ — 2A)a;)/4 > exp(—A„a;)/5} is 
included in An- Since the last inequality and £2 both contain |A„ — A| < r/2 for 7/2 > 0 
small enough, we get 

P(H;)) < P(An > 770) + P(|An - A| > 7/2) + P(A„( 77 .) + 77 ,"^ > 7/3) . 


Clearly the two hrst probabilities in the RHS are 0(77, ^). For n large enough, the last 
probability is less than P(A„( 77 .) > 1^3/2), which is 0(77.“^) by applying Proposition l 5 .lL say 
with r = 77 “^. We thus get El for 77. large enough. 

Proof of bound (j64|l . By (El, and (El, V 3 is dehned as the inverse Fourier 
transform of 


V*{u)=IC{hnU) 


dxCll + dx(l2 Il,n + l2,\ 


(X, W) , 
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where dxtti is a shorthand notation for dai/dx. Using that 


dxdi + 8x0,2 Il,n + /l,' 


< 


. 1 = 1,2 


+ 


dxOi + 8 x 02 


\o„ — o\ 


we obtain that, on the set dehned above, for all i/ G M, 

8 xOi + 8 x 02 




^l,n + ^2,n + £-n 


def 


def 


where, for i = 1 , 2 , we dehne =' e'^"* 8 xOi — {x, ii^) and En = — o\{x, iv). 

Multiplying by 1^^, taking the expectation and applying the Parseval Theorem yield 

2 

\\K\\l 5^supE[lA„^,y + Ml supE[l^„£„ 


dll^lAn] < C2 


i=l 


^eR 


1/GR 


(65) 


where, by (jHI) and (HU and Parseval’s theorem. 


I \K*{Kiy)l 


8 xOi + 8 x02 


{x, in) du < ||iP*||oo 11^112 • 


By (Hi and (Ei, we have 


8xOi - Ji- 


{x, in) < 8 xOi(x, in) - Ii{x, in; A„, P') 


Ii{x, in; Xn, P') - h{x, in; A„, P'n) 


Using this decomposition in Ei^n, the independence of A„, {X'^,Yl), k = 1,... ,n, the fact 
that Var(l(X' < x) < 2 e^^^ and the bound A„ < on An, we get 


< dT2 < E 


1 a„ 


8xOx{x, in) - ii{x, in; A„, P') 


n 


-1 


Using (I49|l and the mean value theorem for bounding the hrst expectation shows that the 
hrst term is 0(l/n) and thus 

supE[£^^„] < (66) 


i/SR 


By (fTj and m, we have 

8x02 - l2,n {x, in) < 8 x 02 {x, in) - l 2 {x, in; A„, P') 

+ hix, in; A„, P') - l 2 {x, in; A„, P'n) 
Using (Ei and (1^ . by the mean value Theorem, we get 

2 I 

< M3 77,“^ 


supE 


8 x 02 {x, in) - l 2 {x, in; A„, P') 

2 ' 

j/SR 















































NONPARAMETRIC INFERENCE ON INDIRECT MEASUREMENTS 


23 


Using (gZl) in the proof of Lemma EH for all 1 / G M, on the set {A„ < r^o}, 

poo 

i2{x,iiy;Xn,P')-i2ix/w;Xn,P'n) < M4 / g{u) \CP' - CP'nlic + iu/w) du , 


where g is an integrable function only depending on c and A. Inserting the three last 
bounds in the definition of £^2,n; we obtain 

supE[l^„£:| J < Cin~^ . 

Comparing m with (0 and m with m, one easily sees that similar argument applies 
for bounding Sn on the set giving 

supE[lA„i^^] < C^n-^ . 


Inserting and the two last displays into (jnni) shows dni- 


Appendix C. Proof of Proposition 15.11 

We have Mi < E[A'] + E[y'] < 00. Denote by Zj {u!i,Ui), i = 1,2 and define the 
function L^{x,y)) Q-{c+ii^)x-wy.^ 

\L^^{x,y)-L^^{x,y)\<g{x,y)\7.i-7.2\i, (x,|/)eM+ (67) 

where g{x,y) max{x,y) and |zi — Z 2 I 1 "= \uJi — uj 2 \ + |pi — P’21- Note that 

A„(IU)= sup \P'nL,-P'L,\. 

x£[-W,W? 

Let N [hU/r]^, where \x~\ denotes the unique integer in [x,x + 1). Then there exists a 
net {zfc}i<fc< 7 v so that 

N 

|-VF,VF]=C [JC(zk,r), 

k=l 

where C'(zk,'r) {z G : |z — Zfc|i < r}. Using this covering in the above expression of 
A„, we get 

An(IU) < max < sup iPA-^z —-P^-hzl • (68) 

1<A:<N j 

Using (EZl) for bounding each term in the max of we get 

sup jP'nL^ - P'L^l 

zeC(zk,r} 

= sup |P'„(Lz — Lzj.) + P'{Lzk ~ -^z) + {P'nLz^. — P'Lzi^)\ 

zeC(zfe,r-) 

<r(P'ng + P'g) + \(Mn-P')Lz, 


(69) 
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Inserting this bonnd in for proving (1^ . it is now snfficient to bonnd P(r[P'„ + 

P']{9) P ^/2) and P(niaxi<fc<Ar \ {P'n — P')Lzf. \ > e/2). Using that P'g = Mi and Markov’s 
ineqnality, we get, if rMi < e:/4, 


P ( r[P'n + P']g > -) = P ( 

^ - 2/ U/2-rMi 


> 1 < 


rMi 


e/2 — rMi 


< 


ArMi 


(70) 


where, in the last ineqnality, we used that l/{x — 1) < 2/x for all x > 2 with x = 
e/{2rMi). Since the LHS is at most 1 and the RHS is more than one when rMi > e/4, 
this inequality holds in all cases, yielding the hrst term in the RHS of ()24j) . We now 
consider the second term in the RHS of (|^ . Since the {{X/, Y/), k >1} are independent 
and identically distributed and |L ^ | < 1 , by using Hoeffding’s inequality (see e.g. Appendix 
6 of I Van Per Waart and Wellnei ( 199 (tIH . we get, for all z, 

- 2 ' 


P ( 


max 

V 1<A;<A 


{P'n-PL)L. 


> 


< AN exp 


ne 

'le 


The proof is concluded by using that N < {W/r + 1)^. To prove inequality ((2SI), let now 
Ly be dehned as L,y(u;x, A) with u = {u,y)] same calculations can 

be done, yielding, for all positive x and A, 


\L^,{u;x,\) - L^2 (u;x, A)| < l[o,x](M)e 


Hu-x) 


■‘"="1 < (?( U )| Z / 1 - Z / 2 | 


where the function g is here dehned as g{{u, y)) *= \y\. Using that P'g < Mi and \Li^\ < 1, 
inequality (I2ni) stems along the same lines as above. 


Appendix D. Proof of Theorem Id. II 
Denote by the integrated workload at time x, that is: 

n = f W{t)dt, (71) 

Jo 

where {W{x), x > 0} is the workload process given in ((T)). Recall that {Ra,, x > 0} denotes 
the on-off process equal to 0 in idle periods and equal to 1 in busy periods (see (j2I)). Dehne 
by p(x, y) the probability: 

p(x, y)=¥ (Aa, = 0, Fa; < 2/) • (72) 

In a hrst step, we calculate the Laplace transform Cp of p using the renewal process 
of the idle and busy periods. Note that this renewal process is stationary. Dehne by 
{Rm n>l} the successive time instants of the end of the busy periods and by {An, n > 1} 
the integrated workload at the end of the busy periods, 

n n 

+ A': = (Zt + xy and 

k=l k=l 

where we have set Rq 0 and Aq *= 0 . 


n > 1, 


(73) 
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Proposition D.l. Under Assumption \(H-l for any {s,p) G such that Iie{s) > 0 
and Re(p) > 0, 


Cp{s,p) = 


XCP'(s,p) 1 

X -^+ 


s + X — XCP'{s,p) p{s + X) p{s + X) 


Proof. The proof is based on classical renewal argnments and the fact that for all integer 
k, the idle period is distribnted according to an exponential distribntion with scale 
parameter A, £x. Note that the event {Sx = 0, < y} may be decomposed as 

{Sx = o,Yx<y} 

= U < T'} U ( y ir; + .Y' < I < < v 




k=l 


(74) 


= {x < T(} U ( IJ {Rn <X<Rn + Zn+1, < I/} j , 

Vn>l / 

where An and Rn are dehned in dZSD. Since Zn+i is independent of these variables, we get 

r+cxD 

p{x, y) — / P(a: — u < Rn < x , An < y) X e“^“ du . 

Writing 

r+oo 

/ P(a; — u < Rn < X , An < y) X e“'^“ du 

Jo 

r+oo 

= P(-Rn < x,An < y) - X / P(i?„ < U - x,An < du , 

Jo 

the proof follows from the identity 


poo /*oo 


lP(-Rn Y X, An < y)e *^e ^'^dxdy = 


'0 JO 


X 


sp \ s + A 


CP'is,p) 


□ 


We will now derive another expression for Cp, nsing standard properties of the Poisson 
process. 

Proposition D.2. Under Assumption \(H-rjjll(H-2)\ for any (s, p) G such that Re(s) > 0 
and Re(p) > 0, 


Cp{s,p) = 


P + OO 


^-(s+A)a; 


p{s + A) p 


exp ( A / e k{x, dv) ) — 1 


dx 
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Proof. Denote by {Mt, t > 0} the counting process associated to the homogeneous Poisson 
process {T^, /c > 0} of the arrivals, more explicitly A/i = l{Tn < t}. By conditioning 

the event {S^ = 0,Yx < y} on the event {A4 = n}, 


p{x, y) = e ^ P(A/; 

n>l 


^{Tj + Xi < X 


k=l 



(75) 


The conditional distribution of the arrival times (Ti,... ,T„) given {A4 = n} is equal to 
the distribution of the order statistics of n independent and identically distributed uniform 
random variables on [0,a;]; hence, for any n-tuple (xi,... ,x„) of positive real numbers, 

P(Di < Xi, . . . , T„ < Xn I A4 = n) = P(f7(l) < Xi, . . . , U(n) < Xn) , (76) 

where {Uk\^=i are independent and identically distributed random variables uniformly 
distributed on [0,x] and U(i) < ■ ■ ■ < f/(„) are the order statistics. Therefore, d7K|) and dZSD 
imply that 


A 


= pl{Ti + X,<x}^^„'^Yk< 

\ k=l 




/ „ n f ^ 1 ^ 

yk<y> ]^D(cixfc 

k=l I k=l ) k=l 


, uykjuuk , 


since the latter integral is invariant by permuting the indexes. An application of the Fubini 
theorem leads to 



where k{x, dy) is the probability kernel defined by 


a(x, dy) 


def 



u)l{u < x}P{du,dy) . 


( 77 ) 


We obtain, for any p such that Re(p) > 0, 


1 A" 


p{x,v)e dv = -h - 


—Xx 


p p ^' n\ 

^ ^ n>l 


k{x, dv)e 


uo 


g—Ax 


P P 


and hence 


Cp{s,p) = 


^H-cxd 


p{s + A) pJo 


^-su^-Xn 


exp ( A J k(x, dv)e J — 1 


exp { ^ J 6 dv) ) — 1 


du . 


□ 
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The proof of Theorem 1,‘f.ll is then a direct consequence of Proposition ID. H and Proposition 
inm and the fact that 

a{x,p) = exp I A / e~^^K{x,dv) 


The result is extrapolated on the line Re(p) = 0 by continuity in p at hxed s such that 
Re(s) > 0. 
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